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The longitudinal viscosity r\i is obtained for a two-dimensional (2D) liquid using a Green-Kubo 
method with a molecular dynamics simulation. The interparticle potential used has the Debye- 
Hiickel or Yukawa form, which models a 2D dusty plasma. The longitudinal r\i and shear rj a viscosities 
are found to have values that match very closely, with only negligible differences for the entire range 
of temperatures that is considered. The bulk viscosity rjb is determined to be either negligibly small 
or not a meaningful transport coefficient, for a 2D Yukawa liquid. 

PACS numbers: 52.27.Lw, 52.27.Gr, 66.20.-d, 83.85.Jn 



I. INTRODUCTION 

The longitudinal viscosity r\i is a transport coefficient 
of interest for fluids pQ. It is the counterpart to the 
better-known transverse viscosity [2] , which is more com- 
monly called the shear viscosity rj s . The latter charac- 
terizes the momentum flux perpendicular to a velocity 
gradient. These viscosities are theoretically predicted to 
be related by [li] 

n d- 1 

Vl = 2— -j— r] a + (1) 

where d is the dimensionality of the system, and r\b is the 
bulk viscosity. The bulk viscosity r/b [5] is also called the 
volume viscosity or expansive viscosity; it is a parameter 
for liquids as well as molecular gases [Bj. Both the shear 
and bulk viscosities appear in the Navier-Stokes equa- 
tion [7] of a fluid. However, as compared with the shear 
viscosity 77,,, the longitudinal viscosity r/i and the bulk 
viscosity 775 are studied less often. 

Physically, these kinds of viscosity characterize energy 
dissipation in a fluid. Bulk viscosity is for energy dissi- 
pation due to compression and rarefaction of a fluid, for 
example in shock waves and high-frequency sound waves. 
Shear viscosity, on the other hand, is for energy dissipa- 
tion due to a gradient in the flow velocity. In the latter 
case, the energy dissipation rate is proportional to both 
the shear viscosity and the square of the velocity gra- 
dient [3 [9] . In the case of a periodic density perturba- 
tion, the energy dissipation rate is proportional to both 
the bulk viscosity and the square of the rate of density 
change [HUE]. 

Unlike shear viscosity, longitudinal and bulk viscosities 
are in general difficult to measure experimentally 1 1 2] . 
This is because the hydrodynamic effects of bulk vis- 
cosity are significant only for rapid time variations, un- 
like shear viscosity which affects flows in easily observed 
ways, even under steady conditions. Ultrasound attenua- 
tion has been described as the only experimental method 
available for measuring bulk viscosity of a fluid [12] . In 
this method, one can obtain the bulk viscosity after sub- 
tracting the ultrasound attenuation contributions from 



thermal conduction and shear viscosity [SJ [T31 H3] ■ This 
method of measuring bulk viscosity is so difficult that 
it has been used only for water and a handful of exotic 
liquids [5], and even for these substances the results for 
the bulk viscosity have large uncertainties. In contrast to 
this difficult experimental situation, however, longitudi- 
nal viscosity and bulk viscosity are mentioned more often 
in the theoretical literature, where it has been calculated 
for example using the Green-Kubo relation J2] [T31l28j . 
using the hydrodynamic limit [29] , or derived using a 
Chapman-Enskog approach 30 . In this paper, we will 
make use of the Green-Kubo approach, which we present 
in Sec. II. We will use the Green-Kubo method to obtain 
r\i and r] s , and for comparison we will use Eq. |l]) to study 
Vb- 

Dusty plasma |31H35| is partially ionized gas contain- 
ing micron-size solid particles, also called dust particles. 
These dust particles are highly charged negatively within 
the plasma by absorbing more electrons than ions, since 
negatively charged electrons have a higher temperature 
than positively charged ions. Due to the shielding pro- 
vided by free electrons and ions in the plasma, the inter- 
action between dust particles in a plasma can be mod- 
eled using a Yukawa or Debye-Hiickel potential [36] , sim- 
ilar to charged particles in a colloidal suspension [37) . 
Because of the high particle charge, dust particles in 
plasmas are strongly coupled (i.e., the potential energy 
between neighboring particles is larger than its kinetic 
energy), so that the collection of dust particles exhibits 
properties of liquids or solids. In laboratory experiments, 
dust particles can be in two-dimensional (2D) or three- 
dimensional (3D) suspensions, depending on the experi- 
mental conditions. In 2D experiments, all dust particles 
are confined in a horizontal plane, with negligible out- 
of-plane motion due to strong confining potentials in the 
vertical direction. The dust particles are immersed in 
a rarefied gas, which applies a much weaker friction to 
moving dust particles, as compared with the case of a 
colloid. The size of dust particles allows imaging them 
directly and tracking their motion, so that various trans- 
port mechanisms can be studied experimentally at the 
particle level. Transport mechanisms that have been 
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studied for dusty plasmas include diffusion [35] , shear vis- 
cosity 39{, and thermal conduction [40 . For 2D dusty 
plasmas, viscosity is generally attributed to dust particle 
scattering arising from interparticle interactions, while 
scattering due to the molecules of rarefied gas is negligi- 
ble, as explained in [UJ . Shear viscosity has been widely 
studied for dusty plasmas, first in simulations for 3D sys- 
tems [2U O [43] , then later in experiments [39] [41] and 
simulations [HI US] for 2D systems, as well as 3D ex- 
periments j2Hl [37] . The longitudinal viscosity and the 
bulk viscosity have been quantified for classical 3D one- 
component plasmas (OCP) [3[[Tlimj[32] using molecular 
dynamics (MD) simulations and gluon plasmas [48 , how- 
ever, it has until now not been quantified for classical 2D 
dusty plasmas, to the best of our knowledge. 

In this paper, we will report a determination of the 
longitudinal viscosity for a 2D Yukawa liquid. We will 
also report an unusual finding regarding the bulk viscos- 
ity: it is either negligibly small or it is not a meaningful 
transport coefficient, for a 2D Yukawa liquid. 



II. GREEN-KUBO RELATIONS FOR 77; AND rj b 

Green-Kubo relations are often used to calculate var- 
ious transport coefficients, such as diffusion [49], shear 
viscosity [UJ SI] > and thermal conductivity [50] . Green- 
Kubo relations arc for equilibrium conditions; they use 
microscopic random motion of particles to determine 
transport coefficients without any macroscopic gradients. 
The longitudinal and bulk viscosities can also be calcu- 
lated using the Green-Kubo relations [U [T4Tf2"B] using 
similar equations as the shear viscosity. The required 
inputs for calculating longitudinal and bulk viscosities 
include time series of particles' positions, velocities, and 
interparticle forces. 

Now we review the three steps for calculating the lon- 
gitudinal viscosity using the standard Green-Kubo re- 
lation [1 H5H2D1 CI2H25]. These Green-Kubo relations, 
which were originally developed for three dimensions 
(d = 3), are adapted here for two dimensions (d = 2) 
by setting the velocity and coordinate in the z direction 
to be zero. 

First, we calculate a diagonal element of the stress ten- 
sor P xx {t), which is defined as 
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Here i and j indicate different particles which all have 
the same mass m, N is the total number of particles, 



(xi,yi) is the position of particle i, x 
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|r» — r^|, and is the interparticle 



potential energy. The positions and velocities of particles 
111 Eq. (2j vary with time, and this accounts for the time 
dependence of P xx (t). The off-diagonal element of the 



stress tensor, P xy (t), can be used to calculate the shear 
viscosity [UJ [44] . Unlike P xy (t) which fluctuates around 
zero, however, P xx (t) fluctuates around a constant level 
P xx (^) ■ 

Second, we calculate an autocorrelation function for 
the fluctuation of P xx (t) using 



C,(t) - ((P xx (t) - P xx (t))(P xx (0) - P xx (t))). (3) 

Here Ci(t) is the stress autocorrelation function. The 
brackets (■ ■ ■) indicate an average over an equilibrium 
ensemble, which in practice we replace by an average 
over different initial conditions. 

Third, we integrate the stress autocorrelation function 
over time to yield the longitudinal viscosity rji [T8] [51] 



1 



Ak B T 



Ci{t)dt, 



(4) 



where A is the area of the 2D system and T is its tem- 
perature. (For a 3D system, A would be replaced by 
the system volume V.) These equations represent the 
Green-Kubo relation for the longitudinal viscosity in 2D 
systems. To improve statistics, we calculate rji twice, 
using P xx as shown above and also using P yy , and we 
average the resulting values of r/i. 

In addition to the longitudinal viscosity, we can also 
calculate the shear viscosity rj s [UJfGJISD] for 2D systems 



using 



N 



i=l 
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'Ji X Viy ^ 2-^1 



Or 



C s (t) = (P xy (t)P xy (0)), 



, (5) 



(6) 



and 



1 



Ak B T J 



C s {t)dt. 



(7) 



The bulk viscosity [1 [H [TJ [5TJ for 2D systems can be 
calculated similarly, using 

1 



Pit) = ^Px,{t) - P xx {t) + P yy (t) - P yy (t)), (8) 



and 



C b (t) = (P(i)P(O)), 



= Ak«T L Cb{t)dt - 



(9) 



(10) 



We will use the same simulation data as the inputs in 
calculations of r/i and ry s . 

It has been questioned theoretically whether transport 
coefficients are meaningful for 2D liquids. This question 
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has been studied theoretically, starting with a 2D hard 
disk system [52] and then liquids with other interparti- 
cle potentials [50]. A transport coefficient is deemed to 
be not meaningful if the corresponding autocorrelation 
function has a long-time tail that decays as slowly as 
1/t, so that the Green-Kubo integral does not converge. 
For a 2D Yukawa liquid, the validity of transport coeffi- 
cients has been discussed in detail in [HI |44j [50l H31 [54] . 
In Sec. IV, we will present our autocorrelation functions 
and discuss whether they have a long-time tail. 

Equations (2fT0 ) are presented in physical units, al- 
though we will perform simulations using dimensionless 
units. Some of the parameters we will use when making 
quantities dimensionless include the area A of the sim- 
ulated system, the areal number density n, the particle 
mass m, the Wigner-Seitz radius a = (nn)^ 1 ^ 2 , a char- 
acteristic plasma frequency [55] to p d = (Q 2 /2-Keoma 3 ) 1 ^ 2 , 
and the particle kinetic temperature T. Here, Q is the 
particle charge. 



III. SIMULATION METHOD 

To model 2D dusty plasmas, we perform equilibrium 
MD simulations using a binary interparticle interaction 
with a Yukawa potential [35] . We integrate the equation 
of motion mrj = —V faj I0r au particles. This equa- 
tion of motion does not include any friction term or any 
Langcvin heating term. Particles are constrained to move 
only within a single 2D plane. Our simulation includes 
N = 1024 particles in a rectangular box with periodic 
boundary conditions to model an infinite system. The 
Yukawa potential is <j>ij = Q 2 exp(— rij / \d) / {^(-orij), 
where Ad is the screening length. We truncate the 
Yukawa potential at distances beyond a cutoff radius of 
24.76 a; this truncation has been justified in [44 . This 
simulated system is essentially the same as a Yukawa 
OCP, except that we constrain the particle to move only 
on a single plane at z = 0. 

Yukawa systems can be described by two dimensionless 
parameters: the coupling parameter T and the screening 
parameter k. They are defined as T = Q 2 /{4Tre ak B T) 
and k = a/Xjj. One can think of T as an inverse temper- 
ature and k as an inverse indicator of density. 

The input parameters in our simulation include k and 
r. We choose a single value of k — 0.5, which is typical 
for 2D dusty plasma experiments [41] , When k = 0.5, the 
melting point of 2D Yukawa system is T ss 142 [56] . To 
study 2D Yukawa liquids over a large temperature range, 
we choose twelve different values of T varying from 140 
(corresponding to a temperature near the melting point) 
down to 2 (corresponding to a much higher temperature). 
The integration time step is in a range of 0.0037 and 
0.037 uJ~d, depending on the choice of T, as in [33]. For 
each value of T, we perform four runs with different initial 
configurations of particles. 
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FIG. 1: The temperature fluctuation during about one half 
of the 10 6 steps used for data analysis for F = 20 and k, — 0.5. 
Here, T _1 is a dimensionless temperature. 



We use a thermostat only for the initial equilibrium 
of our simulation, and not for the data used to calculate 
r/i and ij s . For each simulation run, we first integrate 
10 5 steps using a Nose-Hoover thermostat to approach 
equilibrium at a desired temperature [44] under steady 
conditions. We then turn off this thermostat to integrate 
another 10 6 steps. Only the data in the latter 10 6 steps 
will be used to calculate the viscosities. We use a suffi- 
ciently small time step so that the energy conservation is 
adequately obeyed during the simulation run, as we have 
verified for our simulation data. We measure T, which 
can differ slightly from the desired temperature, using 
the mean square velocity fluctuation. 

Figure 1 shows the time series of the measured temper- 
ature from one of our simulation runs. The temperature 
fluctuates about a steady level during the 10 6 steps sim- 
ulation interval. The temperature fluctuations are due 
to the finite simulation size. The absence of a general 
upward or downward trend in the temperature as a func- 
tion of time is due to our choice of an adequately small 
integration time step. When we report a value for T, we 
use a temperature that was averaged over the f0 6 steps, 
for a given run. 



IV. RESULTS 
A. Longitudinal and shear viscosities 

Using the particles' positions, velocities, and potentials 
from the simulation, we use Eqs. ^ and |5]) to calculate 
the time series of the stress tensor elements. Examples of 
the results for the stress tensor are shown in Fig. 2. All 
our calculations of r\i and r] s will be based on these time 
series. We note that P xx (t) fluctuates about a non-zero 
value, while P X y(t) fluctuates about zero. The source of 
this fluctuation is the microscopic compression and shear 
motion of particles. We find that fluctuations of P xx (t) 
and P xy (t) have comparable amplitudes and time scales. 
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FIG. 2: (Color online). The fluctuations of the stress ten- 
sor elements during about 5% of the 10 6 steps used for data 
analysis for r = 20 and k = 0.5. 



We calculate the autocorrelation function of the time 
series of P xx {t) using Eq. ([3]), and the result is shown 
in Fig. 3. Also shown is the autocorrelation function 
of P xy (t). In panels (a-b), we see the initial decay fol- 
lowed by a small variation around zero. We will use the 
area under these curves to calculate the viscosities, us- 
ing Eq. Q for rji and Eq. Q for r] s . The integral in 
Eq. Q has a finite upper time limit of infinity, which 
we replace with the first zero crossing [JT] , marked as ij 
in Fig. 3. In panel (c) we show the positive portion of 
these autocorrelation functions on a logarithmic scale so 
that the two curves can be distinguished. The two curves 
are almost identical in the initial decay, which leads us 
to expect that the longitudinal and shear viscosities will 
have almost the same values. 

Indeed, we find that the longitudinal viscosity rji and 
the shear viscosity r\ s have almost the same values, for 
the full range of T that we investigate. This is seen in 
Fig. 4(a), where we present rji and rj s determined by per- 
forming the integrals of the stress autocorrelation func- 
tion in Eqs. Q and Q. We find negligible differences 
between them, as shown in Fig. 4(b). 

We also see a minimum in rji as a function of T. This 
minimum matches the minimum in rj Sl which was previ- 
ously studied and explained as being due to a balance of 
kinetic and potential terms in the shear stress [44], [45] . 

In Fig. 4, the data have scatter in both axes. For the 
horizontal axis, the scatter around each value of T in- 
dicates slight differences in the measured temperature, 
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FIG. 3: (Color online). The autocorrelation functions of 
the stress tensor elements, where Ci for longitudinal viscosity 
and C s for shear viscosity, for F = 20 and k = 0.5. The 
same autocorrelation functions are shown with linear axes in 
(a) and (b), and a log axis in (c). Data are normalized by 
parameters defined in Sec. II and III. The first zero crossings 
in (a) and (b) are marked as ti, which will be used to replace 
the upper limits in the integrals, Eqs. Q and Q. 



which can occur due to the absence of a thermostat, as 
discussed in Sec. III. For the vertical axis, the scatter 
corresponds to the random run-to-run variation for the 
obtained viscosity values, i.e., random errors. In addi- 
tion to these random errors, there is a systematic error 
associated with the choice of the upper integral limit. By 
examining the fluctuation of the Green-Kubo integral at 
long times }57) , we determined that this systematic error 
is smaller than the random errors. 
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FIG. 4: (Color online), (a) The longitudinal and shear vis- 
cosities from evaluating Eq. Q and Eq. respectively, (b) 
The difference between rji and r/ a is negligible, for all values 
of r. Using Eq. ( |f I [ ), this result indicates that r]t is negligibly 
small. The minimum in rj s in (a) is known to arise from a bal- 
ance of kinetic and potential terms in the shear stress [441 145] , 
and here we find that a similar minimum appears in rji. 



B. Bulk viscosity 

The negligible difference between rji and r) s that we find 
in Fig. 4 leads us to determine that the bulk viscosity 775 
is much smaller than either 77/ or r^. This conclusion is 
drawn from Eq. ([I]) for our 2D system, which is 



Vb = m - 7] s - 



(11) 



Previous simulations using the Green-Kubo approach 
to obtain the bulk viscosity were mostly for 3D Lennard- 
Jones interparticle potentials and soft-sphere interparti- 
cle potentials [F0HT81 l20l I23H28] , or similar interparticle 
potentials [19 . In some of those simulations, the shear 
viscosity was also calculated [ForlTfl |2"31 |2"41 |2"8] , and it 
was found that the bulk viscosity differs from the shear 
viscosity, with the difference within one order of magni- 
tude. 

Simulations of 3D plasmas [21 El EB [22] provided 
results for the bulk viscosity for both one-component 
plasmas (OCP) and Yukawa OCP. In these simula- 
tions [21 [TH (2TJ [22] , it was found that the bulk viscosity 



is negligible as compared with the shear viscosity, about 
two orders of magnitude smaller or even more. From this 
aspect, it seems that our results for 2D Yukawa system 
that 7/f, -C r/ s are consistent with those previous simula- 
tions in 3D OCP systems. 

We now examine the autocorrelation functions used in 
calculating rn, rj s , and rib in Eqs. (|4j [7] and 10 1 to de- 
termine whether they exhibit a long-time tail. As dis- 
cussed in Sec. II, if the correlation function decays more 
slowly than 1/t, this long-time tail prevents the conver- 
gence of the Green-Kubo integral so that the correspond- 
ing transport coefficient is deemed to be not meaningful. 
In Fig. 5(a-b) we present the correlation functions Ci for 
longitudinal viscosity and C s for shear viscosity, and we 
find that they do not exhibit a noticeable long-time tail 
before the function becomes noisy. However, in Fig. 5(c) 
the correlation function Cb decays more slowly as can be 
seen by comparing to the line drawn with a slope cor- 
responding to a 1/t scaling. This result suggests that, 
within the uncertainties that are inherent in a finite-size 
simulation [50], r\i and r) s are meaningful, but r)b is not. 
It is interesting that the signal-to-noise ratio for Cb, the 
correlation function of the bulk viscosity, is still compa- 
rable to that of C s and C/, even though its amplitude is 
one or two orders of magnitude smaller. Even if 775 were 
meaningful, it would have a small value because Cb in 
Fig. 5(c) is two orders of magnitude smaller than Ci and 
C s . 

We cannot explain in terms of the macroscopic fluid 
equations why the bulk viscosity is either negligibly small 
or not meaningful for this 2D liquid. However, in terms 
of microscopic motion, we can discuss some of the terms 
of the correlation functions. The correlation function 
for the longitudinal viscosity C\ involves only products 
of P xx with a delayed version of itself, and likewise for 
P yy . The bulk viscosity has a different character because 

Eq. ^ also includes cross terms like (P xx (t)P yy (0)). 
In fact, the correlation function for the bulk viscosity, 
Eq. ((9j), can be written as the sum of two terms: Ci/2 
and a cross correlation involving P xx and 
2D Yukawa liquid these two terms almost cancel 



Pyy. For our 



V. SUMMARY 

Molecular-dynamics simulations of a 2D Yukawa liquid 
demonstrate that the longitudinal viscosity is almost the 
same as the shear viscosity, over a wide range of temper- 
ature. These results were obtained using Green-Kubo in- 
tegrals of the appropriate autocorrelation functions. The 
very close match of the values for r]i and r/ s would pre- 
dict, using Eq. 0, that rjb is negligibly small or even 
zero. Indeed, rib might not even be a meaningful trans- 
port coefficient for the system studied here because we 
found that its autocorrelation function exhibits a long- 
time tail, so that the corresponding Green-Kubo integral 
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FIG. 5: (Color online). Correlation functions for (a) rji, (b) n a , and (c) r/b, for F = 140 and k, = 0.5. Here, data are shown 
with log-log axes to allow an identification of any possible long-time tail. We find that only the correlation function Ct, for the 
bulk viscosity in (c) has a significant long-time tail, as seen by a decay that is slower than 1/t. Fhese results, for Y = 140, are 
representative of the other values of Y studied as well. 



diverges. This divergence does not occur for rj s and T)i, as 
they do not have a long-time tail, as judged by our sim- 
ulation. We note that our results are based on a finite- 
size simulation; future larger simulations might be able 
to provide noise-free correlation-function data for longer 
times to further test these conclusions. 

This work was supported by NSF and NASA. 
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